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Abstract In the present study we investigate magnetic reconnection in twisted 
magnetic fluxtubes with different initial configurations. In all considered cases, 
energy release is triggered by the ideal kink instability, which is itself the result 
of applying footpoint rotation to an initially potential field. The main goal of this 
work is to establish the influence of the field topology and various thermodynamic 
effects on the energy release process. Specifically, we investigate convergence 
of the magnetic field at the loop footpoints, atmospheric stratification, as well 
as thermal conduction. In all cases, the application of vortical driving at the 
footpoints of an initally potential field leads to an internal kink instability. With 
the exception of the curved loop with high footpoint convergence, the global 
geometry of the loop change little during the simulation. Footpoint convergence, 
curvature and atmospheric structure clearly influences the rapidity with which 
a loop achieves instability as well as the size of the subsequent energy release. 
Footpoint convergence has a stabilising influence and thus the loop requires more 
energy for instability, which means that the subsequent relaxation has a larger 
heating effect. Large-scale curvature has the opposite result: less energy is needed 
for instability and so the amount of energy released from the field is reduced. In¬ 
troducing a stratified atmosphere gives rise to decaying wave phenomena during 
the driving phase, and also results in a loop that is less stable. 
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1. Introduction 


The solar corona is thought to be heated to temperatures of millions of degrees 
Kelvin by dissipation of stored magnetic energy. The details of the processes 
of energy storage and dissipation remain contentious, and it is likely that the 
corona is heated by a combination of mechanisms (Parnell and De Moortel, 
2012). A very plausible scenario, especially for heating in Active Regions, is 
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that the corona is heated by the combined effect of many small flare-like events 
known as ’’nanoflares” (Parker, 19881. Thus, understanding of flares, especially 
smaller events, contributes to solving the coronal heating problem. It is essential 
to understand the heating of loops, which are the main building blocks of the 
coronal magnetic field (Reale, 20141. 

Twisted magnetic fields are ubiquitous in the solar corona, and twist is asso¬ 
ciated with free magnetic energy. New magnetic flux ropes emerging from below 
the solar surface should already be twisted, and further twisting is produced 
by photospheric footpoint motions with vorticity. In essence, the magnetic fields 
that permeate the solar corona acquire free energy due to the convective motions 
that take place in and around the loop footpoints: i.e., where the field intersects 
the photosphere. There is an increasing body of observational evidence for solar 
flares occurring in twisted coronal loops {e.g. Srivastava et al., 2010 Kuridze 


et al., 2013 Kumar and Cho, 20141, and high resolution observations from HiC 
show untwisting of ’’braided” fields associated with energy release (Cirtain et 
al., 20131. Here, we focus on the energy release within a single coronal loop 
containing magnetic field which is twisted by vortical photospheric motions. 
The primary emphasis is on modelling microflares and similar events, but the 
combined effect of many such events with different magnitudes also provides 


an effective coronal heating mechanism (Browning and Van der Linden, 2003 
Bareford et al. 


2010 , 20111 . 


Previously, numerical simulations have shown that magnetic energy release, 
sufficient for coronal heating above active regions, can occur as a consequence of 
an ideal instability (Browning et al., 2008 Hood et al, 2009[ Botha et al., 2011 


Bareford et al., 20131. As the free magnetic energy increases, the coronal loop 
becomes more and more susceptible to the kink instability (Hood, 19921. This 
instability, although in itself an ideal process, triggers the formation of multiple 
current sheets and leads to magnetic reconnection at many sites within the loop 
volume. The magnetic reconnection causes energy to be released from the field 
and heats the coronal plasma. 

The earlier models of this process start with a field configuration that is 
known to be already unstable (ie.slightly beyond the threshold for linear ideal 
kink instability). In addition, the field is usually modelled as a simple straight 
cylinder where the field strength depends on the distance from the loop axis 
only. In reality, coronal magnetic fields must expand from localised photospheric 
sources - despite the much-discussed phenomenon of “constant cross-section”. 


which now seems more likely to reflect the plasma emission (Klimchuk, 2000 


Peter and Bingert, 20121. Furthermore, coronal loops are also inevitably curved. 

These simplifications may affect the dynamics and energetics of twisted loops, 
artificially restricting or exaggerating the amount of magnetic energy released. 
As well as the field geometry, features involving atmospheric physics, which 
are expected to influence how well magnetic energy is thermalised within the 
loop volume, are also usually ignored. For example, most previous work uses a 
simplified adiabatic energy equation, although the effects of thermal conduction 
in cylindrical unstable loops have been considered by Botha et al. (2011). 

Recently, Gordovskyy et al. (2014) investigated the kink instability and mag¬ 
netic reconnection in more realistic configurations with large-scale curvature 
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(ie.loop-like fluxtubes) and atmospheric stratification. Although the key focus 
of that study is particle acceleration, it has revealed some features not present 
in idealised cylindrical models, such as the loop contraction, asymmetric current 
distribution and variation of the energy release along the loop. Hence, the loop 
topology and various thermodynamic effects can substantially affect the energy 
release process and, therefore, they need to be investigated in more detail. 

The aim of this paper is to consider heating in twisted coronal loops, involving 
more realistic models, and, along the way, we explore the effects of removing var¬ 
ious idealisations and simplifications which have been adopted in previous works. 
Does the basic mechanism of an ideal kink instability triggering the formation of 
fragmented current sheets, leading to significant dissipation of magnetic energy, 
persist in a realistic loop model? If so, how are the dynamics and energetics 
affected by the magnetic field geometry and the plasma physics incorporated 
within the modelling? Therefore, we investigate a family of four loop models 
which allow us to compare and contrast the effects of field geometry, both loop 
curvature and field expansion/convergence. In contrast to previous models which 
assume an initial unstable equilibrium, the twisted fields are established by 
rotating the footpoints of the field. Furthermore, in the case of the more realistic 
curved loop models, we explore the effects of incorporating thermal conduction, 
and also embed the loop in a gravitationally-stratified atmosphere. In order to 
study the effect of stratification, we introduce an additional model, with the 
atmosphere stratified in both density and temperature, thus representing a loop 
whose lower layers are embedded in the upper chromosphere. 

The paper is structured in the following manner. Section describes the dif¬ 
ferent loop configurations used for the nonlinear magnetohydrodynamic (MHD) 
simulations. Section [^presents the numerical code and Section [^focuses on how 
results change as the loop model is made more and more realistic: e.g., how the 
energy release is affected by the addition of curvature, or how the results change 
when a stratified atmosphere is used. Finally, in the last section, the results are 
summarised and our conclusions are given. 


2. Loop configurations 

The primary goal of this study is to investigate the effect of field geometry on the 
energy release in twisted loops. Hence, our simulations cover four separate field 
configurations, see Figure [l] In models A and B, we consider straight fluxtubes 
using a domain with dimensions x = ±5 Mm, y = ±5 Mm, z = ±10 Mm, with 
the fluxtubes initially orientated along the z-axis. The first fluxtube (A) has 
a constant cross section, whereas model B has a field converging towards the 
footpoints: the field strength at the footpoints is twice the value at the centre 
(or apex) of the fluxtube. The purpose of loop B then is to show the influence 
of field convergence with respect to the kink instability and subsequent heating. 
Loops C and D have more realistic configurations, featuring large-scale curvature 
and are orientated differently compared to the straight loops. The curved loop 
domains have dimensions x = ±10 Mm, y = ±10 Mm, z = 0 Mm - 20 Mm, with 
both footpoints residing at the z = 0 boundary, representing the chromosphere. 
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Figure 1. Initial magnetic field geometry for loops A (top left), B (top right), C 
(bottom left) and D (bottom right). 
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Figure 2. Left, the footpoint driving function for O.lftw ^ f ^ O.Qtt-w Right, the Alfven 
speed along the central field line of loops C (green) and D (blue) — 6^ = 90 is the 
position at the loop apex. 


Loop C is constructed such that the field converges towards the footpoints in 
a manner comparable to loop B (both have Bftp/Bapx = 2): the differences in 
results between these two loops should therefore reveal the impact of curvature. 
The footpoint positions for loop D are the same as those for loop C, and so 
these two loops exhibit a similar level of curvature; however, loop D is given five 
times the level of footpoint convergence (Sftp /Bapx = 10): the intention here is 
to continue to explore the impact of field convergence, but within the context of 
loop curvature. 

The initial field for loop A is uniform and has a z-component only, while in 
loops B, C, and D the initial field is constructed using two point sources located 
outside the model domain: 


B{r) 




r- S 2 \ 
\r - S2|3y ’ 


( 1 ) 


where r is a position vector within the simulation volume, Si and S 2 are the 
positive and negative point sources; Bg is a simple scaling constant and is given a 
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specific value such that the axial field strength at the footpoints is some multiple 
of the field strength at the apex. In model B these point sources are located 
outside the ^min and Zmax boundaries, while in C and D they are both located 
below Zmin boundary. Table gives the point source coordinates and scaling 
factors required to implement the converged field geometries. 


Table 1. The point source coordinates and scaling terms required 
for field convergence (see eqn.Q. 


Loop 

Si 

S2 

Bs 

Bftp 

Bapx 

A 

n/a 

n/a 

n/a 

1 

1 

B 

(0, 0, -21.4) 

(0, 0, 21.4) 

114 

1 

0.5 

C 

(0, -7.5, -27.1) 

(0, 7.5, -27.1) 

1480 

1 

0.5 

D 

(0, -7.5, -4.2) 

(0, 7.5, -4.2) 

17.9 

1 

0.1 


The density and temperature are initially uniform in models A, B, C and D: 
n = 1.2 X 10^® m“^ and T = 4 x 10^ K. In addition, model D is considered with a 


gravitationally stratified atmosphere (as model D*), see Section 4.4 


For all loops discussed in this paper, the magnetic field is under the line-tied 
condition (ry = 0, i.e.dB/dt = 0 when u = 0) at the “foot-point” boundaries. 
Furthermore, all loop simulations begin with a potential field. The magnetic 
twist necessary for instability is created by rotating the plasma located at the 
footpoints; rotation vortices are, for the curved loops, centered at a; = 0, y = ±7.5, 
z = 0. The rotational driving varies in space and time (Figure]^ left): 


ve{r,t) = rujo 


1 - tank (^) 


X tanh ( 

^sw 


1 — tanh 


{ht-) 


( 2 ) 


where r is the radial distance from the vortex centre, wq = 0.015 is a scaling 
factor, Tfr = 0.6 Rq is the footpoint radius, rfb = 0.05 Rq is the footpoint boundary 
thickness, tsw = 20 Ia is the switching time, and ftw is the characteristic twisting 
time. There is a slightly different arrangement for the straight loops: the vortex 
centres are located at (0,0, ±10) and the two footpoints are driven in opposite 
directions. Note, the dimensionalising factors are Ro = l Mm and tA = 0.35 s (see 
Section for further details). The important parameters are Wq and the 
product of these two terms, multiplied by two, gives the total twist after footpoint 
rotation has stopped: e.i?., for ftw = 750 the twist is ~67r. Each loop is twisted 
for a time sufficient to cause an instability. The Poynting flux generated by the 
driving adds energy to the field; if the driving is continued through the unstable 
phase it becomes difficult to assess the energy released as the loop relaxes to a 
lower energy state. Thus, ftw changes according to the loop: the values of the 
twisting time for the five loops are 750 (A), 650 (B), 400 (C) and 600 (D). 

The level of field convergence determines how the Alfven speed varies with 
height. For loop A there is no convergence and va ~ 2800 km s“^ for all z, which 
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is also true for the footpoints of loops B-D. Hence, the peak driving velocity 
(~ 21kms“^) is sub-Alfvenic (Figure]^ right) for all four loops. 


3. Numerical code 


All numerical simulations were performed using LARE3D - 3D Lagrangian 


remap MHD code (Arber et al., 20011. This code is based on a Lagrangian 
remap scheme: the Lagrangian part, which is done using a second-order accurate 
predictor-corrector method, deforms the grid such that it moves with the plasma. 
The advantage of this technique is that additional physics, such as thermal 
conduction and shock capturing, can easily be incorporated into the code. The 
remap stage involves the mapping of the plasma properties {e.g., density, velocity, 
magnetic field) back to the original Cartesian grid; monotonicity is preserved 


through the use of Van Leer (1997) gradient limiters. 


LARE3D solves the following resistive MHD equations. 


d 


^ = 
dt 

pv) = 


-V-(pu), (3) 

—V • ( pvv ) -f — ( V X b\ X B — VP — pgz + V ■ a (4) 

Un V / 


dB 


= V X 


(»xB) - Vx 


— (pe) = -V-(peu) - PV v + -|- V q 

with specific energy density, 


+ ea, 


(7 - 1) P ’ 


( 5 ) 

( 6 ) 

( 7 ) 


and current density 


J = —V xB. 
Mo 


( 8 ) 


where p is the mass density, v is the plasma velocity, B the magnetic field, P the 
thermal pressure, t] is the resistivity, 7 = 5/3 is the ratio of specific heats, po is the 
magnetic permeability, and q is the conductive heat flux. Radiation is ignored in 
this study, but it is not expected to have a significant impact, since the radiative 
timescale is much longer than the magnetic relaxation timescale. The gravity is 
ignored {g = 0) in models with uniform atmosphere (A-D), but included in model 
D* with atmospheric stratification. The last terms of equations and feature 
tensors, and are required to simulate shock heating; these terms are defined in 
Section 13.21 
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We normalise the variables in the MHD equations using reference values 
suitable for a coronal active region: 

r* p* ^ B* 

tro Po ^0 

where asterisks denote the unnormalised MHD variables, Rq = 1 Mm, pq = 
2 X 10“^^ kg m“^ and Bq = 4.47 x 10“^ T. Other variables are expressed as 
follows; 


L - ^ 


t* V* „ P* 

tA VA Pq 


where va = BqI y/poPo is the reference Alfven speed, tA = Ro/ va is the reference 
Alfven time, and Pq = Bq /pq is the reference pressure. The specific energy den¬ 
sity, current density and resistivity (e, J and 77 ) also have reference variables 
that can be expressed in terms of i?o, Po and Bq: 


Co = 


3. 

PoPo 


= vi, 


Jo = 


Bn 


PqRq ’ 


770 = poRqva- 


Hence, for the chosen values of Rq, po and Hg, ~ 0.36 s, va ~ 2800 kms“^ and 
770 « I.Itt X 10® n m. 

In all four cases, the simulation features two stages (see also Gordovskyy 
et al., 20131: twisting until the fluxtube becomes kink-unstable and magnetic 


relaxation after the kink instability. The first stage is done using ideal MHD: 
Ohmic dissipation and conduction are absent from the energy equation. The 
footpoint driving eventually induces a kink instability that initiates a release of 
magnetic energy. Just before this point, the simulations are restarted with the 
resistive and conductive terms now included in the MHD equations. Any heating 
will reduce density and thereby increase the local Alfven speed, which then leads 
to a decrease in the time step under the CFL condition. For the converged loops 
the magnetic field is strongest at the footpoints, and so the driving will inevitably 
create the highest currents at these locations too. The burst of Ohmic heating 
caused by the switch on of resistive MHD can be so intense as to cause cavitation 
at the z boundaries; the result is that the time step becomes too small for the 
simulation to make progress. Fortunately, this issue can be avoided through the 
use of thermal conduction (Section |3.1[ ), assuming that the plasma below the 
footpoints is maintained at a constant cold temperature equivalent to 4000 K. 

Throughout the loop volume there is a background resistivity of 77 b = 3 x 10“^, 
except when the current reaches or exceeds a threshold (jcrit = 2 ), at which point 
an anomalous resistivity, rjc = 0.002, is applied. The lower the value of jcrit the 
greater the percentage of the loop interior that will be assigned an anomalous 
resistivity when the simulation is restarted in resistive mode. We chose a current 
threshold of two, since this results in 10% of loop A contributing to Ohmic 
heating immediately after the restart. 

The computational domain is a 3D staggered grid: physical variables are not 
calculated at the same place for each cell in the domain, which improves numer¬ 
ical stability and allows conservation laws to be included in the computation. 
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There are some differences between the straight and curved loops as regards grid 
dimensions and limits. The straight loop simulations are run at a grid resolution 
of 128^ X 256, whereas for the curved loops it is 256^ x 512. The differences in 
grid volume and resolution between the straight and curved loop simulations 
mean that, along the x and y axes, the straight loops are better resolved by a 
factor of two; however, the curved loops have double the resolution along the z 
axis. 

3.1. Thermal conduction 

The conductive heat flux (equation is implemented using Braginskii parallel 
conduction (Braginskii, |1965[ ) which becomes isotropic where the magnetic field 
is weak: 


9 


B 


\B\ 


VT ]B 


■ 

min 


\B\ 


■ VT, 


( 9 ) 


where k. = kq (kq = 10“^^) and 6min = 0.001. Conduction results in the trans-| 
fer of internal energy across the grid and is calculated implicitly over half the 
time step, using successive over relaxation. Before equation|^can be applied the 
code variables must be converted from normalised units to SI units. This is done 
by amending the kappa constant. 


Kg 


1 

LopI pI 



( 10 ) 


Then the changes in internal energy must be converted back into normalised 
units before the code can proceed. 


3.2. Shock resolution 


LARE3D uses shock viscosity (Wilkins, 19801 to capture the heating effect of 
shocks, see the last term of the energy equation Q. It is expressed as the product 
of the rate of strain tensor, 


1 I dvi dvj 


^ii ... 


2 \ di 


+ 


d: 


and the shock tensor. 


^ij — P ^ms T ^2 ^ ^ ’ 


( 11 ) 


( 12 ) 


where Cms = \/cg + is the magnetosonic speed, 4 is the distance across a grid 
cell in the direction normal to the shock front, s is a similarly localised strain 
rate (the subscripts i and j denote the different spatial coordinates), and the 
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shock viscosity coefficients 1 ^ 1 = 0.1 and 1^2 = 0.5 are constants (real viscosity is 
set to zero). 

The method of using shock viscosity to represent heating due to shocks is 
known to yield physically valid results. Recently, Bareford and Hood (20151 
conducted a detailed investigation of shock handling within LARE3D. They have 
shown that, for kink-unstable loop simulations, LARE3D viscous heating is con¬ 
sistent with Petschek reconnection and slow-mode shocks. Unlike the anomalous 
resistivity, controlled by the externally-defined critical current, the dissipation 
due to shock viscosity is defined internally, by the velocity field in the simulation 
domain. 


4. Results of numerical simulations 


This section presents the results taken from the loop simulations (A-D). Some 
properties are given in SI units; these are length, temperature and density. 
Velocities are expressed in either kms“^ or Mms“^. Times are given in units of 
the Alfven time, whereas energies and currents are given in normalised units. 
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Figure 3. Left, the magnetic twist in units of tt as a function of radial distance from 
footpoint centre for the straight loops, A (black) and B (red). Right, the same plots 
but for the curved loops, C (green) and D (blue). The twist values were determined 
numerically just before instability onset, see time labels. The twis t thres hold for an 
ideal kink instability (0«2.497r) as calculated by Hood and Priest (19791 is given by 
the dashed horizontal line. 
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Figure 4. The change in magnetic energy normalised by the initial value, £'b(0) (left) 
and the natural logarithm of the kinetic energy (right) for loop A (black) and for loop 
B (red). Resistive MHD was switched on at t = 750tA for loop A and t = 650tA for 
loop B , i.e., just before instability. Energies are volume intergrated. 
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Figure 5. The cumulative Ohmic and shock heating (left) and the change in internal 
energy (right) for loop A (black) and for loop B (red). Energies are volume intergrated 
and normalised by the initial magnetic energy, i?b(0). 
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Figure 6. Left, the percentage of the grid assigned anomalous resistivity for loop 
A (black) and for loop B (red). Right, the percentage of the cross sectional area 
(0 < r < 1) assigned anomalous resistivity at t = 750 tA. 
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Figure 7. The current density magnitude (dark regions correspond to the highest 
currents) over the x-y plane at t = 751tA and « = 2.4 (loop A, left) and 2 = 0.6 (loop 
B, right). 
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Figure 8. The temperature over the y-z 
(middle) and 1200 tA (bottom) for loop A 
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plane (a; = 0) at times t = 750 (top), 760 
left) and for loop B (right). 
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Figure 9. The average velocity parallel to the z-axis at three times for loops A (left) and 
B (right). The velocity component Vz is averaged over 0 < r < 0.8 for each z coordinate. 


4.1. Kink instability 

First, we check that the driving phase has generated sufficient twist consistent 
with the ideal kink instability. Figure is created by following field lines from 
specific points on the postive footpoint. The starting points are distinguished by 
the distance from the footpoint centre, located at (0,0,-10) for the straight loops 
and (0,-7.5,0) for the curved; and the magnetic twist (LBg/rBz) is calculated 
numerically by keeping track of how many times a field line wraps around the 
initial loop axis before it reaches the negative footpoint. The twist measurements 
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were taken at a time close to instability onset. It can be clearly seen that the 
level of twist for both straight loops exceeds 2 . 497 r, the result obtained by Hood 
and Priest (1979) which is applicable to a straight loop of aspect ratio ten 
{ = L/R): hence, this value should only be treated as a necessary condition for 
instability. At radii less than 0.4, the magnetic twist is approximately b-Tir: this is 
consistent with the average twist on the instability threshold (where the twist is 
single-signed) reported by Bareford et al. (2011). However, the curved loops are 
significantly less twisted: in particular, the twist profile for loop D is never greater 
than Stt and is only above the Hood and Priest threshold when 0.35<r<0.4. 
These results suggest that curvature exerts a destabilising influence: i.e., the 
average absolute twist at instability onset is lower than it is for the straight 
loops. 


4.2. Footpoint convergence 

In order to investigate the importance of footpoint convergence we compare the 
results for loops A and B. The former begins with a uniform straight held of 
I S I = 1 , whereas the latter begins with a held that has less energy, since | B \ 
only rises to one at the footpoints. 

Figures]^ and [^present the changes in volume-integrated energies and, with 
the exception of the kinetic energy logarithm, these terms have been normalised 
by the initial magnetic energy, i?b(0), in order to allow a meaningful comparison 
between the two straight loops. The vortical footpoint driving (equation]^ leads 
to instability for both cases; however, the converged straight loop (B) is the hrst 
to achieve instability. For this reason, resistive MHD is switched on earlier, see 
the vertical dash lines in Figure]^ 

Proportionally, the driving phase adds more azimuthal held to loop B, even 
though it is driven for less time: furthermore, after t = 650tA, loop B undergoes 
two bursts of energy release. The hrst occurs, as expected, at the end of the 
driving phase, then the loop stabilises for 300 fA before undergoing a slightly 
greater burst of energy release (similarly, loop A shows a quasi-stable period 
centred on 800 fA)- By the end of the simulation, loop B has released more than 
double the magnetic energy released by loop A. As a result, the converged loop 
produces higher levels of kinetic and internal energy. Incidentally, the energy 
released by loop B increases only marginally (~ 8 %) should the driving phase 
be extended to t = 750 Ia (as is the case with loop A): the instability is delayed 
until this later time after which the magnetic energy declines to a relaxed state, 
with no intervening stable period. 

It is expected that the driving should generate a linear increase in azimuthal 
field over time: hence, the increase in magnetic energy should show a depen¬ 
dence. This relationship is conhrmed for loop A by the dashed line in Figure]^ 
(left), although, there is clearly a faster increase in magnetic energy when the 
field converges at the footpoints. After the driving phase, the natural logarithm 
of the kinetic energy rises again when resistive MHD is switched on (the gradient 
of Inifkin is twice that of 7 , the growth rate of the instability, since the velocity 
of the perturbed plasma is proportional to e'*'*). The wave-like form present from 
the start of the simulations is a consequence of the fact that the driving speed is 
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supersonic; however, by t = A00tA, these oscillations have diffused away, which 
indicates that the field is now evolving through a sequence of equilibria (Mellor 
et al., 2005). Overall, the numerical dissipation, measured as a percentage of the 
total initial energy, is 0.11% for both straight loops. 

The levels of Ohmic and shock heating are roughly equal for the unconverged 
loop (A), but when the field is converged the shock heating becomes steadily 
greater as the simulation progresses: in fact, Ohmic heating only increases during 
the two periods of magnetic energy decline. 

Figurej^shows the percentage of the grid assigned anomalous resistivity (left), 
which, of course, is zero until the start of the resistive phase. The plot on the 
right gives the location of supercritical currents as a function of z: the percentage 
value refers to the portion of a circle (r = 1) whose origin follows the z-axis that 
has 7^ = 0.002 (Figure]^ left). 

Inevitably, the driving creates high currents at the footpoints — these are not 
shown in Figure]^ (right) in order to reveal how rj changes around the loop apex. 
Nevertheless, for loop A, 90% of the loop volume assigned anomalous resistivity 
occurs within —9 < z < 9, whereas loop B has 82% of the anomalous resistivity 
within this range. 

We now investigate the current structures existing shortly after instability 
onset (Figures]^ and 10). Similar to previous simulations of kink-unstable loops 
(Browning et al, |2008 Hood et al., 2009 Botha et al., 2011 Bareford et al., 
2013[), the current sheets in the cylindrical models A and B form a helical ribbon 


structure around the kinking fluxtubes. In addition, there is a strong currect 
along the fluxtube axis (x = y = 0). This current structure has elliptic cross- 
section with its main axis rotating along the z-axis; the total rotation angle just 
after the kink is approximately the same as the total magnetic twist angle. In 
the model A (with initially cylindrical fluxtube) the current density is neary 
uniform along z-axis: it is slightly higher in the central region of the fluxtube 
(ie.around z = 0), while in the fluxtube with converging field (model B) it also 
has high current densities near foot-points (see also Gordovskyy and Browning, 


2011 20121 . 


Next, we examine how heat is distributed within the two straight loops at 
three times, t = 750tA (instability onset), 760(immediately after the instabil¬ 
ity) and 1200<A (when the field has relaxed to a lower energy state), see Figure 
1^ Leading up to the instability (top row), hot plasma forms shells around the 
fluxtubes, apparently, corresponding to the current density concentrations. At 
this stage the temperatures are also comparable. However, immediately after the 
instability, temperatures for loop A increase by a factor of four, with the highest 
temperatures forming along the axis. The radial spread of heating at certain 
locations along the z-axis reveals where the loop has kinked. Contrastingly, loop 
B shows the strongest heating near the footpoints, as well as in the shell around 
the fluxtube; the temperatures are roughly double those before the instability. 
As the loops relax, thermal conduction acts to smooth out the temperatures, 
yielding an average value of 2 MK. 

Finally in this section, we look at how Vz varies along the loop axis (Figure 
1^. The velocity is averaged over the loop cross-section (0 < r < 0.8) for the same 
times as the temperature plots. Before the instability (solid line), the stongest 
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axial flows for loop A are away from the apex. These flows diminish as the 
instability progresses, resulting in a residual left to right flow by the time loop 
A has relaxed. Loop B has a more complicated flow pattern, consistent with 
heating sources located at 2 = ± 7.5. Over time the flow becomes chaotic, but 
the end result is similar to loop A, albeit with a higher rightward flow. 


4.3. Loop curvature 


The obvious difference between curved loops (models C and D) and cylindrical 
fluxtubes is the geometry of current density distribution. In the straight loop 
(model A) currents are nearly uniformly distributed along the loop, while in the 
converging fluxtube (model B) the current density is, as one would expect, higher 
near the footpoints. However, in both models the structure has high degree of 
cylindrical symmetry. At the same time, in curved loops (models C and D) the 
current distributions do not have cylindrical symmetry. Thus, just before the 
instability occurs, the current is concentrated in a thin shield above the loop top 
(see Figure [l0|. In addition, there are current concentrations close to footpoints 
due to the field convergence. 

Loop D has the most highly converged footpoint field, which means the field 
strength at the apex is the weakest (all loops have the same footpoint field 
strength). The result of this is that loop D has by far the smallest volume- 
integrated field energy — over ten times smaller than the total field energy 
for the other curved loop (C). Hence, the energy plots for these two models are 
strikingly different. Figuresandare not normalised by the initial integrated 
field energy, since other plots would appear almost flat by comparison (Figure 
right). 

As with the straight loops, increasing the convergence, increases the rate at 
which the loop accrues free magnetic energy during the driving phase; but, within 
the context of curvature, greater convergence appears to delay the instability by 
some 200 tA- Loop D has to be driven for at least 7001 a before an instability 
will occur, whereas loop C encounters an instability at 400 ^a- 

In addition, the two curved loop simulations appear to differ as regards the 
nature of the instability. Loop C is consistent with the kink-like instabilities seen 
for the straight loops: there is a swift drop in magnetic energy coincident with 
rises in heating (both Ohmic and shock) and internal energy. Furthermore, the 
unstable phase is also accompanied by peaks in kinetic energy. Loop D on the 


other hand, shows a more gradual drop in magnetic energy (Figure 121, which 


does not correspond to the rise in Ohmic heating; instead the similarly gradual 
increase in internal energy is caused by shock heating. Surprisingly, the kinetic 
energy remains high even as loop D settles into a lower energy state. Driving 
loop C for the same amount of time as loop D (be., ttw ~ 6001 a) does little 
to alter how the magnetic energy changes during the simulation. The instability 
occurs at roughly the same time as before (400 ^a), but the energy released is 
halved: the continued driving interferes with how the instability plays out. 

To understand the impact of curvature, we compare the normalised change in 
magnetic energy from the straight loops with loop C (Figure [T^ left). Comparing 
the red and green plots, we see that curvature has throttled the rate at which 
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driving adds magnetic energy to the loop. At the same time however, large-scale 
curvature has made the loop more susceptible to instability: a loss of magnetic 
energy occurs much earlier. 

In the right panel of Figure we add the result from loop D to show 
how continuing to increase footpoint convergence almost completely cancels the 
effect of curvature. The fall in magnetic energy still happens earlier than for the 
straight loops, but the growth in magnetic energy before instability far outstrips 
that seen in the other simulations. The normalised increase in magnetic energy 
is eight times higher than that for loop B. 

We now turn our attention to showing whether or not the results from loop 
D are indeed consistent with a kink instability. Although, Figure]^ (right) shows 
that loop D is sufficiently twisted (stronger than the loop C), we also see that the 
changes in magnetic energy and the kinetic energy logarithm (Figure |11[ blue 
lines) are different from the forms seen for other kink-unstable loops. Specifi¬ 
cally, the decline in magnetic energy is noticeably slower and not accompanied 
by a rapid change in kinetic energy. It seems that increasing the convergence 
(compared to loop C) has prolonged the transition to a low-energy state. 

However, the magnetic fieldline plots in Figures [l4| and 15 (right columns) do 
agree with the form expected for a kink instability: initially, the fieldlines are 
tightly twisted around the loop axis, and then, as the instability proceeds, the 
fieldlines untwist. Loop C undergoes an internal kink instability, since the apex 
height (zwS.S) remains almost constant throughout the simulation, and the 
heating (a mixture of Ohmic and shock) is concentrated around the apex during 
the unstable phase. As the loop relaxes, conduction combined with continued 
heating results in a near-uniform temperature of around 1.5 MK. The tempera¬ 
ture plots for loop D do suggest an increase in apex height. On the other hand, 
the corresponding fieldline plots indicate that the apex rise is temporary and is 
no longer evident once the loop has relaxed. Figurej^also shows that the (mostly 
shock) heating is concentrated along the legs of the loop, and, for the scaling 
used here {Bq = 44.7G), results in sub-MK temperatures. There is some Ohmic 
heating at instability onset (t = 600tA), but it is confined to the footpoints. 

Comparison of all four models shows that the ratio of viscous to Ohmic 
dissipation changes from ^ 1 for a non-converging fluxtube (see also Hood et 
ai, 2009) to as much as ~ 10 in the strongly converging loop in model D. While 
the contribution of the Ohmic and viscous heating in models A and B is nearly 
uniform along the fluxtubes, in curved fluxtubes in models C and D the Ohmic 
heating is prevalent very close to footpoints (where the current density is high), 
while in the middle of the loops, where velocitites and velocity gradients are 
high, the magnetic energy dissipates through the “viscosity” channel. Hence, 
the enhanced viscous dissipation in converging fluxtubes is most likely caused 
by shocks formed due to low magnetic field and, hence, low Alfven velocities 
near the fluxtube centres. 

Strong shock dissipation in model D can explain the delay of the kink instabil¬ 
ity. Since the shock viscosity effectively acts as an additional Ohmic dissipation 
(see Section 3.21, the rates of helicity injection and magnetic energy accumulation 
(due to the footpoint rotation) are lower. However, as far as the reconnection 
stage is concerned, slower energy release is due to the field convergence: lower 
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magnetic field near the looptop and, hence, lower current densities, reduce the 
Ohmic heating rate, which is proportional to . 


4.4. Stratified atmosphere 


Loop C 
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Figure 10. Current density distributions jnst before the kink instability in models C 
and D. Upper panels show loop mid-planes (x = 0), lower panels show central cross- 
sections (j/ = 0). Current density scales are shown in units of jo = 3.6 x 10“^ A m~^. 


The loops mentioned so far all have a density and temperature that are ini¬ 
tially uniform, where n = 1.2 x 10^® m“^ and T = 4 x 10^ K. In order to investigate 
the effect of a stratified atmosphere, we also have loop D*, identical to loop D, 
except that it has a non-uniform atmosphere based on the work by Gordovskyy 
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Figure 11. The change in magnetic energy (left) and the natural logarithm of the 
kinetic energy (right) for loop C (green) and for loop D (blue). Resistive MHD was 
switched on at t = 400 Ia for loop C and t = 600 Ia for loop D, i.e., just before instability. 
Energies are volume intergrated. 
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Figure 12. The cumulative Ohmic and shock heating (left) and the change in internal 
energy (right) for loop C (green) and for loop D (blue). Energies are volume intergrated. 



Figure 13. Left, the volume-integrated change in magnetic energy normalised by the 
initial value, Eb(0) for loops A (black), B (red) and C (green). Right, the same plot 
but with the change in magnetic energy for loop D (blue). The weightings used for loops 
A and B also take into account the difference in grid volume between the straight and 
curved loop simulations: the latter being sixteen times the former. 


et al. (2013): 


p{z) = Pi exp 


-(zsh + z) 


Zl 


+ p 2 exp 


-(Zsh + z) 


Z2 


(13) 


where pi = 3.34 x 10 ® kg m ^ is the photospheric density, p 2 = 2 x 10 kg m ^ 
is the mean density in the corona, zi=0.25Mm is the density scale height 
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Figure 14. Loop C, the temperature over a loop cross section (left column) taken 
at three times, instability onset (top), during the instability (middle) and the relaxed 
state (bottom). The configuration of the field lines, alongwith the pattern of Ohmic 
(orange) and shock (purple) heating (right column), are shown for the same times. 


between the photosphere and the transition region, Z 2 = 50 Mm is the coronal 
scale height and 2 :sh = l-5Mm simply allows the density profile to be shifted 
horizontally. 

The initial atmosphere is shown in Figure [16] — the temperature profile fol¬ 
lows from hydrostatic balance, and so a gravitational term (p g{z)) is added to 
the force equation ([^. Between z = 0 and z = 4 the temperature increases by 
more than two orders of magnitude (from around 4000 K to 0.83 MK), while the 
particle number density drops from 10^° to 10^^ m“^. This region represents the 
chromosphere and transition region. The rest of the domain {z > 4) represents 
the corona, where the temperature and density are more or less constant. It is 
appropriate to add an atmosphere to loop D, since the increase in thermal pres¬ 
sure towards the footpoints is matched by the hundredfold increase in magnetic 
pressure. 

Increasing the density reduces the Alfven speed, which has the unwelcome 
consequence of violating the constraint for direct current heating. At the foot- 
points, the driving speed (Figure left) is now twice the Alfven speed, and 
so the perturbations generated by the driving occur too frequently for the loop 
to have time to settle into an equilibrium. This issue is easily rectified if the 
driving factor, wq, is reduced by an order of magnitude. Hence, ua(- 2 = 0) is now 
five times the driving speed, which compares well with observed values for active 
region field strengths (0.1 T), photospheric flow speeds (lkms“^) and densities 
(4 X 10-‘^kgm-3). 

The stratified atmosphere used for loop D* is not in equilibrium in the pres¬ 
ence of thermal conduction; however, such an equilibirum can be achieved by 
allowing the simulation to evolve without driving for 2000 <a- Only then do we 
begin the driving phase with ttw = 4300set such that the driving ramps down 
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Figure 15. Loop D, the format of this figure follows that used for Figure |14[ 
At t = 600tA, the temperature plot is saturated: near the footpoints, the highest 
temperatures are 6MK. 



Figure 16. The initial density (black) and temperature (grey) with height {z) for loop 
D*. 


at 6300^A- Note, no restart is required after the driving phase, since in order 
for the atmosphere to acquire its initial equilibirum, the simulation must have 
resistive MHD and conduction switched on from the start. 

The following energy plots begin from when the driving is started. 
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Figure 17. The change in magnetic energy (top left), the kinetic energy (top right), 
the cumulative Ohmic and shock heating (bottom left), and the change in internal 
energy (bottom right) for loop D*gravitation. Energies are volume intergrated. 
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Figure 18. The magnetic twist in units of tt as a function of radial distance from 
footpoint centre for loops, D (blue) and D* (green). The format of this figure follows 
that used for Figure]^ 


Figure 17 (top left) shows that the driving takes 4300 tA to induce an in¬ 
stability. Comparison with loop D (Figure [TT|), reveals that the presence of an 
atmosphere requires less buildup in magnetic energy for instability onset, and 
so the resulting energy release is half that seen for loop D. This also means 
that adding a stratified atmosphere reduces the twist required for the instability 
onset: just before the kink instability loop D* has a lower twist for r < 0.5 
compared to loop D. The kinetic energy plot for loop D* exhibits the expected 
peak at the time the instability occurs, but there is also a wave pattern that, 
during the driving phase, decays in amplitude. Figure 17 (bottom left) shows 
that shock heating continues to dominate over Ohmic heating. Internal energy 
undergoes a continual increase, the gradient of which steepens as the loop goes 
unstable. The heating before the instability is possibly connected to the decaying 
wave form seen in the kinetic energy plot. Although disguised by the resolution 
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of the time axis, the relaxation phase is even more extended than it is for loop 
D. The relaxation time is roughly 1000tA, which equates to a dimensionalised 
time of 355 s if one uses the scalings given in Section 


5. Summary and conclusions 


We have shown that the application of vortical driving to the footpoints of an 
initially potential field does lead to a kink instability, for straight and curved 
loops, with or without a stratified atmosphere. The rate at which the driving 
adds free magnetic energy to the loop is certainly enhanced when the field is 
made to converge at the foopoints. Thus, although converging loop B requires 
more energy to become unstable compared to its cylindrical counterpart (A), 
the buildup of energy is sufficiently fast that, for the same driving profile, loop 
B achieves instability around 100 <a before loop A does. When comparing the 
curved loops (C and D), we again see that the strongly converging loop D has 
proportionally more free energy when instability occurs. Although, unlike the 
pair of loops A and B, the strongly converging loop D requires higher twist to 
become unstable compared to the weekly converging loop C. However, this case 
might be not representative, as the substantial viscous heating could affect the 
loop evolution prior to the instability. 

The strongly converging loop D* (which is similar to the loop D, but embed¬ 
ded into the gravitationally stratified atmosphere) reaches the kink instability 
with a total twist nearly equal that in the loop C. Furthermore, taking into 
account the difference in twisting speed, it appears that they require the same 
amount of foot-point rotation to get the kink instability. Hence, one can conclude 
(by comparing models A and B) that the field convergence reduces the stability 
of a twisted loop in cylindrically symmetric magnetic fluxtubes. However, the 
convergence does not seem to influence the stability of the fluxtubes with large- 
scale curvature. These conclusions with regard to the loop stability should be 
used with caution, as the stability can be influenced by many effects, which can’t 
be examined using only five test models. 

The large-scale curvature affects the geometry of magnetic reconnection. The 
current distribution loses its cylindrical symmetry, with strong currents formed 
in a thin shield above the loop top. Although, physically it is represented by two 
effects - twist reduction and fluxtube expansion (due to the reconnection with 
ambient field) - the reconnection is more localised around the central part of the 


loop, leading to topological evolution descibed by (Gordovskyy et al, 20141. 

As regards heating, loops A and C show equal amounts of Ohmic and shock 
heating, whereas loops B, D and D* are mostly heated through shocks. The shock 


heating term (Section 3.21 in the energy equation becomes more significant when 
higher energies are required to drive the loop unstable. The dominance of shock 
heating is also true when there is a stratified atmosphere. The photospheric 
density for D* is some five orders of magnitude higher than the initially uniform 
density used for loop D. For this reason, the driving speed is reduced by a factor 
of ten, so that it is five times lower than the Alfven speed at z = 0, consistent 
with observed values. These two speeds are sufficiently close for wave phenomena 
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to be visible in the kinetic energy profile. The decaying nature of these waves 
suggests that some heating of the plasma is occurring during the driving phase 
(Figuretop right), which perhaps destablises the loop. 

The important methodological implication from our result is that the shock 
viscosity needs to be taken into account as an additional Ohmic dissipation. 
This issue is particularly important for numerical experiments with localised 
resistivity effects, as shocks would affect effective spatial resistivity distribution. 
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